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We present approximate analytical solutions to the Hamil- 
tonian and momentum constraint equations, corresponding 
to systems composed of two black holes with arbitrary lin- 
ear and angular momentum. The analytical nature of these 
initial data solutions makes them easier to implement in nu- 
merical evolutions than the traditional numerical approach of 
solving the elliptic equations derived from the Einstein con- 
straints. Although in general the problem of setting up initial 
conditions for black hole binary simulations is complicated by 
the presence of singularities, we show that the methods pre- 
sented in this work provide initial data with h and loo norms 
of violation of the constraint equations falling below those 
of the truncation error (residual error due to discretization) 
present in finite difference codes for the range of grid resolu- 
tions currently used. Thus, these data sets are suitable for 
use in evolution codes. Detailed results are presented for the 
case of a head-on collision of two equal-mass M black holes 
with specific angular momentum 0.5Af at an initial separa- 
tion of lOM. A straightforward superposition method yields 
data adequate for resolutions oi h = M/4, and an "attenu- 
ated" superposition yields data usable to resolutions at least 
as fine as h — M/8. In addition, the attenuated approximate 
data may be more tractable in a full (computational) exact 
solution to the initial value problem. 

PACS number(s): 04.70.Bw,04.25.Dm 



I. INTRODUCTION 

The computation of gravitational wave production 
from the interaction and merger of compact astrophysi- 
cal objects is an analytical and computational challenge. 
These calculations would be able to provide both a pre- 
dictive and an analytical resource for the gravitational- 
wave interferometric detectors such as LIGO Virgo 
and GEO600 |^ soon to be online. We concentrate on 
the case of binary black hole mergers Q . 

Besides ours, the other major efforts to numerically 
simulate black hole coalescences are being carried out by 
the AEI-WashU-NCSA g and the Cornell-Illinois col- 
laborations Q. The resolutions currently under consid- 
eration by the AEI-WashU-NCSA effort are as fine as 
h — M/5 in a computational domain of with 385"^ mesh- 
points, where h denotes the grid spacing. Our compu- 



tations are carried out at a similarly modest resolution 
h — Af/4, on 161'^ domains. The goal of all these groups 
is to perform eventually simulations in domains of up 
to 1000"^ grid-points, which will allow finer resolution 
and/or large physical domain sizes. 

Regardless of resolution, in order to carry out such sim- 
ulations, we must construct initial data sets representing 
binary black hole systems. These data sets should not 
only satisfy the Einstein constraints but also carry the 
desired physical content. Binary black hole initial data 
sets are the subject of this paper. 

Under the 3-1-1 formulation of General Relativity, the 
construction of initial data requires solving the Hamilto- 
nian and momentum constraints: 

R+^K^- A/ a; = (1) 



and 



V,A/ - -V,K - 



(2) 



respectively. Above, R is the 3-dimensional Ricci scalar 
constructed from the spatial 3-metric gij, and 
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K, 



A, 



(3) 



is the extrinsic curvature tensor. K and Aij are the trace 
and trace- free parts of Kij, respectively. Covariant dif- 
ferentiation with respect to gij is denoted by V^. Space- 
time indices will be denoted by Greek letters and spatial 
indices by Latin letters. 

The systematic solution of the constraint equations is 
due to Lichnerowicz, Choquet-Bruhat, York and others 
Q . This procedure involves freely specifying a metric gij , 
a traceless extrinsic curvature Aij and the trace K. It 
then introduces a conformal factor 0, and a scaling rule 
to determine the physical and Kij : 



9ij 



Ai-j — 
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as functions of the coordinates (defined both in the phys- 
ical and in the trial conformal space). Solution consists 
of simultaneously determining (Eq. (^ becomes an el- 
liptic equation for (p) and correcting the longitudinal part 
of Aij . Typically, one introduces a vector potential to 
accomplish this latter task and the problem consists of 
four coupled elliptic equations for (j) and 0|. 
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II. KERR-SCHILD SLICING 



Our work is based on descriptions of black holes in 
ingoing Eddington-Finkelstein (iEF) coordinates. This 
choice is motivated by the fact that, in these coordinates, 
surfaces of constant time "penetrate" the event horizon. 
Foliations that penetrate the horizon facilitate the exci- 
sion of the singularity from the computational domain. 
The essence of black hole excision is the removal of the 
singularity while preserving the integrity of the space- 
time accessible to observers outside the black hole. As 
originally suggested by Unruh this is only possible 
if the excised region is fully contained within the event 
horizon, thus the need to have access to the interior of 
the black holes. 

Other work, such as that of Briigmann Q for the 
generic 3-dimensional code, does not use excision, and 
instead solves K — Q initial data in a conformally flat 
background Here, we keep K non-zero to main- 

tain close similarities to the analytically known single 
Kerr-Schild black hole. 

The 4-dimensional form of the Kerr-Schild spacetime 

ds^ = -dt^ -f dx^ + dy^ + dz^ + 2H{x°'){hdx^f (4) 

describes isolated single black holes. Here the scalar 
function H has a known form, and Ix is an ingoing null 
vector congruence associated with the solution. For in- 
stance for the Schwarzschild solution, H{x") = Al/r, and 
l\ = (l;a;*/r), i.e. an inward pointing null vector with 
unit spatial part. From the fact that Eq. (Q) is an ex- 
act solution to the Einstein equations, one can write the 
3-metric, and the momenta Kij associated with Eq. (jj). 

Initial data setting for multiple black hole spacetimes 
using the method described by Matzner, Huq, and Shoe- 
maker (MHS) 10] begins by specifying a conformal spatial 
metric which is a straightforward superposition of two 
Kerr-Schild single hole (spatial) metrics: 

(jijdx^dx^ = Sijdx'dx^ + 2 iH{x'^){iljdx^Y 

+2 2Hix''){2ljdx' f . (5) 

The fields marked with the pre-index 1 (2) correspond to 
an isolated black hole with specific angular momentum 
ai (a2) and boosted with velocity Vi (V2). The super- 
position of the conformal momenta is defined as follows: 
The extrinsic curvature for a single hole (say hole 1) 

iK,j = {id,P, + id,(3j - 2 iV'l^iPk - l^tg^j)/{2 la) , 

is converted to a mixed-index object. 

The trace of K is calculated as the sum of the corre- 
sponding traces: 



and the transverse-traceless part of the extrinsic curva- 
ture as 

Ai^,K^+2K^-\5^k . (6) 

MHS propose the use the metric and extrinsic curva- 
ture so defined as a conformal metric and extrinsic cur- 
vature to solve the coupled elliptic system [ p^ . 

III. APPROXIMATE SOLUTIONS 

Here we take a different approach. After all, when two 
black holes are widely separated we expect almost-linear 
superposition to hold, so the data setting as specified 
above, without the elliptic solution for (p and w', should 
lead to only small errors in the widely separated case. 
We show here that even for interestingly close separation 
scenarios, the superposition errors are small, and in fact 
both the ^1 and loo norms are smaller than those of the 
truncation error (the discretization error in the calcula- 
tion) in simulations at currently accessible computational 
resolutions. In all cases here we present the results for 
head-on collisions (which are simpler to display) but very 
similar results are found at similar separations for non 
head-on data. 

Figure |l| shows plots of the norms of the residuals of 
the Hamiltonian and momentum constraints for a head- 
on collision using initial data {g^, Aij} provided by Eq. 
(^) and Eq. (^ , respectively. The parameters of the data 
are Mi = M2 = M = 1, |ai| = |a2| = a = 0.5M along 
the z axis, and the holes are boosted against each other 
in the x direction with velocity v = 0.5. The residuals 
are defined as the absolute value of 

H = R+'^k^-A^'A^' , (7) 

Af, = V,i/ - ^V.X ; (8) 

where again "~" denotes analytic quantities evaluated us- 
ing the approximate solutions (^) and (^). We see from 
figure ^ that the violation of the constraints is primar- 
ily confined to a region near each hole. This is due to 
the fact that Eq. (|^) is the sum of terms which scale 
as 0(Af^/r^ x Al/d), where r is the coordinate distance 
to the singularity and d is the coordinate separation be- 
tween holes. For an isolated black hole, these terms can- 
cel each other out exactly, since the Kerr-Schild metric 
is an exact solution of the constraint equations. How- 
ever, that is not the case for the spacetime metric and 
extrinsic curvatures provided by Eqs. (^) and (|^). The 
presence of the perturbation of the second hole in the 
vicinity of the first hole destroys the balance between 
these terms, causing the right hand side of Eq. (0) to 
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scale as 0{M'^ /r'^ x M/d). As in our computational ap- 
proach to evolving the spacetime, we mask the interior 
of the black hole. In this example, we excise the points 
inside a sphere of radius a + h centered at the hole, where 
the absolute value of the specific angular momentum a is 
also the radius of the Kerr ring-like singularity and h is 
the grid spacing (for figure |l|, /i = M/4) . Hence, figure |l| 
plots precisely the error on the computational domain. 

Every finite difference code will have a corresponding 
truncation error associated with discretization. In most 
of the finite difference codes under development, this er- 
ror scales as 0{h?), second order accuracy; therefore, 
given a second order finite difference discretization of the 
constraints, the truncation errors are obtained from 

Htr^R+^K^ -A^^A^' (9) 

M,tr = VjA^' - '^V.K . (10) 

Above, "~" is used to denote finite difference discretiza- 
tion, and it is understood that these finite difference ex- 
pressions are to be evaluated using the exact, not ap- 
proximate, solutions of the constraints. Of course, these 
exact solutions are not yet available ||l^. In order to 
circumvent this problem, we obtain an estimate of the 
truncation errors from the absolute value of 

H = R+\k^-A^'A^'-H (11) 

m = Vj A/' - ^y,K - M, (12) 

instead. Although these truncation errors correspond to 
a modified form of the constraints, these errors should not 
be substantially different since the structure of the terms 
involving finite difference operators was not changed, 
only analytic functions {H and Mi) have been added. 
Notice that for an isolated hole, the Kerr-Schild spatial 
metric and extrinsic curvature satisfy exactly the con- 
straint equations [H = Mi = 0), leaving us with an er- 
ror exclusively related to the truncation error. Figure |^ 
shows the residual H (diamonds) and truncation error H 
(circles) in the Hamiltonian constraint in a region around 
the rightmost hole of figure ^. The truncation error H 
was obtained in a mesh with grid spacing h = Af/4. 

In order to compare the global degree of satisfaction 
of the initial value equations by the methods described 
here, we calculate the lea (maximum) and li (average of 
the absolute value) norms of the Hamiltonian and mo- 
mentum constraints. We do this over the axis joining 
the black holes from x — —lOM to x — lOM, excluding 
the points inside two segments of length 0.75Af, centered 
around each hole, as indicated above. As shown in Ta- 
ble I, lines 1 and 4, both the li and the l^o norms of 
the superposition error are less than those of the M/4 
truncation error. However figure ^ makes it clear that 
the truncation error is not uniformly (pointwise) larger 



than the superposition error. We find in practice that at 
M/4:, superposition yields adequate behavior in the evo- 
lution. Figure |^ shows a similar diagram to figure ^ for 
the momentum constraint. In this case the superposition 
method can be used to provide data with error pointwise 
uniformly smaller than the truncation error. 

IV. ATTENUATED SUPERPOSITION METHOD 

In order to reduce the residual errors, we propose a 
variation of the superposition method that preserves the 
simplicity of being analytical. Essentially, the method 
consists of multiplying "attenuation" functions into the 
recipe of the previous section. The new approximate met- 
ric (jij, trace K and tensor A/ take the form: 

gfj = 5y + 2 iB{x'') ii?(x'=) ih ilj 

if ' = ,B li^/ + 2B 2AV' - \5.,'k^ . 

The purpose of the attenuation function B is to mini- 
mize the effects due to a given hole on the neighborhood 
of the other hole. For instance, the attenuation function 
iB is unity everywhere except in the vicinity of hole- 
2 where it rapidly vanishes so the metric and extrinsic 
curvature there are effectively that of a single black hole. 
The attenuation functions with this property can be con- 
structed in a number of different ways. An example of 
an attenuation function is 

iB^l- e-^^l"^ , (13) 

where 




P = \J 2l'^{x - 2xY + (y - 2yY + {z ~ 2zY ■ (15) 



Here 22^*, a and 27 denote the coordinate location, spe- 
cific angular momentum and boost factor of hole-2, re- 
spectively. <T represents a free parameter of the attenu- 
ation function. An expression similar to (|l3|) for 2-B is 
obtained by reversing the labels. 

The attenuation function ( p^ was chosen amongst a 
few different types for its simplicity and better perfor- 
mance. Since the constraints involve second order deriva- 
tives, it is important to pick attenuation functions that 
vanish up to second order derivatives, as in (^3|), so 
"pure" single black hole solutions can be obtained in the 
neighborhood of each hole. 
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V. RESULTS 

Figures || and ^ (|| and H) show the Hamihonian (mo- 
mentum) constraint for a region around the rightmost 
hole of figure [|. Figure ^ shows a comparison between 
the truncation error for h — A//4 (empty circles) and the 
superposition initial data (full diamonds), while figure^ 
shows the same for attenuated initial data (full squares). 

First we note in figure |^ that in some areas (near 
X — 2M, for instance) the violation of the constraints for 
these attenuated data is greater than that correspond- 
ing to truncation error. However, even in these areas 
the violation is smooth and small in absolute value, as 
opposed to the unbounded behavior shown by the su- 
perposition data near the singularity, and the divergent 
behavior present at the singularity has disappeared with 
the use of attenuated data. This is due to the fact that, 
having cancelled the influence of the presence of the sec- 
ond hole, the fields and become the fields corre- 
sponding to an isolated black hole at the location of the 
singularity. 

The results are presented in Table |[ The momentum 
constraint h norm represents the average of the three 
components of Eq. (|^) and the loo the maximum value 
of them. These values are compared to the norms of the 
truncation error in the constraints caused by the finite 
difference stencils for grid spacings h — M/4, M/8, and 
M/16. These results show that with these norms, the 
superposition initial data violate the Hamiltonian and 
momentum constraints below the truncation error corre- 
sponding to grid spacings as fine as /i = M/4, while the 
attenuated data extend this range to at least h = M/8. 
Note also that, for both the superposition and the atten- 
uated method, the momentum constraint equations are 
satisfied well within the truncation error for grid spacings 
as fine as h — M/16. Finally Table | (rows 2, 3) shows 
the approximate second order convergence between reso- 
lutions h = M/8 and h = M/16. The convergence rate is 
approximately second order because even at h = M/16, 
finite h effects affect significantly the convergence. 

VI. CONCLUSIONS 

The initial value problem for black hole binary systems 
is currently approached by solving numerically a set of 
elliptic equations derived from the Hamiltonian and mo- 
mentum constraint. Due to the presence of singularities 
in the fields, these equations are complex and difficult to 
handle. 

In this article we presented analytical methods that 
can provide initial value data for black hole binary sys- 
tems. While these solutions are only approximate, they 
satisfy the Hamiltonian and momentum constraint equa- 
tions to the level of accuracy (evaluated with the li and 



Zoo norms) of the truncation error present in finite differ- 
ence codes for a given range of grid resolutions. Because 
attenuated data are so much smoother than the superpo- 
sition, and because the principal difficulties in computa- 
tional simulations arise from sharp gradients and diver- 
gent values, which the attenuation method eliminates, we 
are confident that at least for the particular example of a 
head-on collision of two equal-mass black holes at a sep- 
aration of lOM, thess data could be used in evolutionary 
codes with grid spacings as fine as ft. = M/8. 

The analytical nature of these methods makes them 
simpler to implement than the numerical approach. Fur- 
thermore, a more comprehensive study of the attenuation 
functions may extended the range of resolutions for which 
this approximation is valid. 

The attenuation method offers the advantage that it 
cancels the divergent behavior of the constraints near the 
central singularity. The suppression of contributions near 
the black holes (individual holes have exactly zero residu- 
als) should simplify the exact problem, since this method 
provides exact inner boundary conditions for the elliptic 
solver. 
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FIG. 2. Hamiltonian constraint near the rightmost hole for 
the system of figure ^ The diamonds show the Hamiltonian 
constraint residuals \H\ for the superposition method. The 
circles are the estimation |-ff | of the truncation error present 
in a second order finite difference code {h = M/4), defined by 
Eq. 



♦ ♦ 

* »IIII H II H IIItllll H ltlll H « * 



FIG. 1. Hamiltonian and momentum constraint residuals 
for a head-on collision of two equal-mass M black holes. The 
vertical axis shows the residuals of the Hamiltonian constraint 
|-ff| and the x component of the momentum constraint |Ma;| 
defined by Eqs. (^^- The black holes are separated lOM 
and boosted toward each other at 0.5c. Both black holes have 
specific angular momentum 0.5M pointing perpendicular to 
the colliding direction. The momentum residual corresponds 
to the component parallel to the colliding axis. 




FIG. 3. Momentum constraint (x component) for the sys- 
tem of figure ^. The diamonds show the momentum con- 
straint residuals \Mx\ for the superposition method. The 
circles are the estimation \Mx \ of the truncation error 
{h = M/A), defined by Eq. (Fa). 
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FIG. 4. Hamiltonian constraint for the system of figure |l| 
The squares present the Hamiltonian constraint residuals \H\ 
for the attenuation method. The circles are the estimation 
\H\ of the truncation error present in a second order finite 
difference code {h = M/4), defined by Eq. (0). 
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FIG. 5. Momentum constraint (x component) for the sys- 
tem of figure ^. The squares show the momentum constraint 
residuals \Mx \ for the attenuation method. The circles are the 
estimation \Mx\ of the truncation error {h = M/4), defined 
byEq. (0). 



TABLE I. Truncation errors (first three rows) and residu- 
als (last two rows) of the Hamiltonian and momentum con- 
straints. 



Method 


Ham. 


Constr. 


Mom. 


Constr. 






h 




h 


h = M/4 


4.240 


0.25f 


20.20 


0.605 


h = M/8 


0.757 


0.022 


3.098 


0.065 


h = M/16 


0.f98 


0.004 


0.742 


0.012 


without attenuation 


1.643 


0.097 


0.235 


0.003 


with attenuation 


0.126 


0.049 


0.041 


0.005 
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